home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dlasy2.f < prev    next >
Text File  |  1996-07-19  |  12KB  |  383 lines
  1.       SUBROUTINE DLASY2( LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR,
  2.      $                   LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO )
  3. *
  4. *  -- LAPACK auxiliary routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     October 31, 1992
  8. *
  9. *     .. Scalar Arguments ..
  10.       LOGICAL            LTRANL, LTRANR
  11.       INTEGER            INFO, ISGN, LDB, LDTL, LDTR, LDX, N1, N2
  12.       DOUBLE PRECISION   SCALE, XNORM
  13. *     ..
  14. *     .. Array Arguments ..
  15.       DOUBLE PRECISION   B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ),
  16.      $                   X( LDX, * )
  17. *     ..
  18. *
  19. *  Purpose
  20. *  =======
  21. *
  22. *  DLASY2 solves for the N1 by N2 matrix X, 1 <= N1,N2 <= 2, in
  23. *
  24. *         op(TL)*X + ISGN*X*op(TR) = SCALE*B,
  25. *
  26. *  where TL is N1 by N1, TR is N2 by N2, B is N1 by N2, and ISGN = 1 or
  27. *  -1.  op(T) = T or T', where T' denotes the transpose of T.
  28. *
  29. *  Arguments
  30. *  =========
  31. *
  32. *  LTRANL  (input) LOGICAL
  33. *          On entry, LTRANL specifies the op(TL):
  34. *             = .FALSE., op(TL) = TL,
  35. *             = .TRUE., op(TL) = TL'.
  36. *
  37. *  LTRANR  (input) LOGICAL
  38. *          On entry, LTRANR specifies the op(TR):
  39. *            = .FALSE., op(TR) = TR,
  40. *            = .TRUE., op(TR) = TR'.
  41. *
  42. *  ISGN    (input) INTEGER
  43. *          On entry, ISGN specifies the sign of the equation
  44. *          as described before. ISGN may only be 1 or -1.
  45. *
  46. *  N1      (input) INTEGER
  47. *          On entry, N1 specifies the order of matrix TL.
  48. *          N1 may only be 0, 1 or 2.
  49. *
  50. *  N2      (input) INTEGER
  51. *          On entry, N2 specifies the order of matrix TR.
  52. *          N2 may only be 0, 1 or 2.
  53. *
  54. *  TL      (input) DOUBLE PRECISION array, dimension (LDTL,2)
  55. *          On entry, TL contains an N1 by N1 matrix.
  56. *
  57. *  LDTL    (input) INTEGER
  58. *          The leading dimension of the matrix TL. LDTL >= max(1,N1).
  59. *
  60. *  TR      (input) DOUBLE PRECISION array, dimension (LDTR,2)
  61. *          On entry, TR contains an N2 by N2 matrix.
  62. *
  63. *  LDTR    (input) INTEGER
  64. *          The leading dimension of the matrix TR. LDTR >= max(1,N2).
  65. *
  66. *  B       (input) DOUBLE PRECISION array, dimension (LDB,2)
  67. *          On entry, the N1 by N2 matrix B contains the right-hand
  68. *          side of the equation.
  69. *
  70. *  LDB     (input) INTEGER
  71. *          The leading dimension of the matrix B. LDB >= max(1,N1).
  72. *
  73. *  SCALE   (output) DOUBLE PRECISION
  74. *          On exit, SCALE contains the scale factor. SCALE is chosen
  75. *          less than or equal to 1 to prevent the solution overflowing.
  76. *
  77. *  X       (output) DOUBLE PRECISION array, dimension (LDX,2)
  78. *          On exit, X contains the N1 by N2 solution.
  79. *
  80. *  LDX     (input) INTEGER
  81. *          The leading dimension of the matrix X. LDX >= max(1,N1).
  82. *
  83. *  XNORM   (output) DOUBLE PRECISION
  84. *          On exit, XNORM is the infinity-norm of the solution.
  85. *
  86. *  INFO    (output) INTEGER
  87. *          On exit, INFO is set to
  88. *             0: successful exit.
  89. *             1: TL and TR have too close eigenvalues, so TL or
  90. *                TR is perturbed to get a nonsingular equation.
  91. *          NOTE: In the interests of speed, this routine does not
  92. *                check the inputs for errors.
  93. *
  94. * =====================================================================
  95. *
  96. *     .. Parameters ..
  97.       DOUBLE PRECISION   ZERO, ONE
  98.       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  99.       DOUBLE PRECISION   TWO, HALF, EIGHT
  100.       PARAMETER          ( TWO = 2.0D+0, HALF = 0.5D+0, EIGHT = 8.0D+0 )
  101. *     ..
  102. *     .. Local Scalars ..
  103.       LOGICAL            BSWAP, XSWAP
  104.       INTEGER            I, IP, IPIV, IPSV, J, JP, JPSV, K
  105.       DOUBLE PRECISION   BET, EPS, GAM, L21, SGN, SMIN, SMLNUM, TAU1,
  106.      $                   TEMP, U11, U12, U22, XMAX
  107. *     ..
  108. *     .. Local Arrays ..
  109.       LOGICAL            BSWPIV( 4 ), XSWPIV( 4 )
  110.       INTEGER            JPIV( 4 ), LOCL21( 4 ), LOCU12( 4 ),
  111.      $                   LOCU22( 4 )
  112.       DOUBLE PRECISION   BTMP( 4 ), T16( 4, 4 ), TMP( 4 ), X2( 2 )
  113. *     ..
  114. *     .. External Functions ..
  115.       INTEGER            IDAMAX
  116.       DOUBLE PRECISION   DLAMCH
  117.       EXTERNAL           IDAMAX, DLAMCH
  118. *     ..
  119. *     .. External Subroutines ..
  120.       EXTERNAL           DCOPY, DSWAP
  121. *     ..
  122. *     .. Intrinsic Functions ..
  123.       INTRINSIC          ABS, MAX
  124. *     ..
  125. *     .. Data statements ..
  126.       DATA               LOCU12 / 3, 4, 1, 2 / , LOCL21 / 2, 1, 4, 3 / ,
  127.      $                   LOCU22 / 4, 3, 2, 1 /
  128.       DATA               XSWPIV / .FALSE., .FALSE., .TRUE., .TRUE. /
  129.       DATA               BSWPIV / .FALSE., .TRUE., .FALSE., .TRUE. /
  130. *     ..
  131. *     .. Executable Statements ..
  132. *
  133. *     Do not check the input parameters for errors
  134. *
  135.       INFO = 0
  136. *
  137. *     Quick return if possible
  138. *
  139.       IF( N1.EQ.0 .OR. N2.EQ.0 )
  140.      $   RETURN
  141. *
  142. *     Set constants to control overflow
  143. *
  144.       EPS = DLAMCH( 'P' )
  145.       SMLNUM = DLAMCH( 'S' ) / EPS
  146.       SGN = ISGN
  147. *
  148.       K = N1 + N1 + N2 - 2
  149.       GO TO ( 10, 20, 30, 50 )K
  150. *
  151. *     1 by 1: TL11*X + SGN*X*TR11 = B11
  152. *
  153.    10 CONTINUE
  154.       TAU1 = TL( 1, 1 ) + SGN*TR( 1, 1 )
  155.       BET = ABS( TAU1 )
  156.       IF( BET.LE.SMLNUM ) THEN
  157.          TAU1 = SMLNUM
  158.          BET = SMLNUM
  159.          INFO = 1
  160.       END IF
  161. *
  162.       SCALE = ONE
  163.       GAM = ABS( B( 1, 1 ) )
  164.       IF( SMLNUM*GAM.GT.BET )
  165.      $   SCALE = ONE / GAM
  166. *
  167.       X( 1, 1 ) = ( B( 1, 1 )*SCALE ) / TAU1
  168.       XNORM = ABS( X( 1, 1 ) )
  169.       RETURN
  170. *
  171. *     1 by 2:
  172. *     TL11*[X11 X12] + ISGN*[X11 X12]*op[TR11 TR12]  = [B11 B12]
  173. *                                       [TR21 TR22]
  174. *
  175.    20 CONTINUE
  176. *
  177.       SMIN = MAX( EPS*MAX( ABS( TL( 1, 1 ) ), ABS( TR( 1, 1 ) ),
  178.      $       ABS( TR( 1, 2 ) ), ABS( TR( 2, 1 ) ), ABS( TR( 2, 2 ) ) ),
  179.      $       SMLNUM )
  180.       TMP( 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
  181.       TMP( 4 ) = TL( 1, 1 ) + SGN*TR( 2, 2 )
  182.       IF( LTRANR ) THEN
  183.          TMP( 2 ) = SGN*TR( 2, 1 )
  184.          TMP( 3 ) = SGN*TR( 1, 2 )
  185.       ELSE
  186.          TMP( 2 ) = SGN*TR( 1, 2 )
  187.          TMP( 3 ) = SGN*TR( 2, 1 )
  188.       END IF
  189.       BTMP( 1 ) = B( 1, 1 )
  190.       BTMP( 2 ) = B( 1, 2 )
  191.       GO TO 40
  192. *
  193. *     2 by 1:
  194. *          op[TL11 TL12]*[X11] + ISGN* [X11]*TR11  = [B11]
  195. *            [TL21 TL22] [X21]         [X21]         [B21]
  196. *
  197.    30 CONTINUE
  198.       SMIN = MAX( EPS*MAX( ABS( TR( 1, 1 ) ), ABS( TL( 1, 1 ) ),
  199.      $       ABS( TL( 1, 2 ) ), ABS( TL( 2, 1 ) ), ABS( TL( 2, 2 ) ) ),
  200.      $       SMLNUM )
  201.       TMP( 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
  202.       TMP( 4 ) = TL( 2, 2 ) + SGN*TR( 1, 1 )
  203.       IF( LTRANL ) THEN
  204.          TMP( 2 ) = TL( 1, 2 )
  205.          TMP( 3 ) = TL( 2, 1 )
  206.       ELSE
  207.          TMP( 2 ) = TL( 2, 1 )
  208.          TMP( 3 ) = TL( 1, 2 )
  209.       END IF
  210.       BTMP( 1 ) = B( 1, 1 )
  211.       BTMP( 2 ) = B( 2, 1 )
  212.    40 CONTINUE
  213. *
  214. *     Solve 2 by 2 system using complete pivoting.
  215. *     Set pivots less than SMIN to SMIN.
  216. *
  217.       IPIV = IDAMAX( 4, TMP, 1 )
  218.       U11 = TMP( IPIV )
  219.       IF( ABS( U11 ).LE.SMIN ) THEN
  220.          INFO = 1
  221.          U11 = SMIN
  222.       END IF
  223.       U12 = TMP( LOCU12( IPIV ) )
  224.       L21 = TMP( LOCL21( IPIV ) ) / U11
  225.       U22 = TMP( LOCU22( IPIV ) ) - U12*L21
  226.       XSWAP = XSWPIV( IPIV )
  227.       BSWAP = BSWPIV( IPIV )
  228.       IF( ABS( U22 ).LE.SMIN ) THEN
  229.          INFO = 1
  230.          U22 = SMIN
  231.       END IF
  232.       IF( BSWAP ) THEN
  233.          TEMP = BTMP( 2 )
  234.          BTMP( 2 ) = BTMP( 1 ) - L21*TEMP
  235.          BTMP( 1 ) = TEMP
  236.       ELSE
  237.          BTMP( 2 ) = BTMP( 2 ) - L21*BTMP( 1 )
  238.       END IF
  239.       SCALE = ONE
  240.       IF( ( TWO*SMLNUM )*ABS( BTMP( 2 ) ).GT.ABS( U22 ) .OR.
  241.      $    ( TWO*SMLNUM )*ABS( BTMP( 1 ) ).GT.ABS( U11 ) ) THEN
  242.          SCALE = HALF / MAX( ABS( BTMP( 1 ) ), ABS( BTMP( 2 ) ) )
  243.          BTMP( 1 ) = BTMP( 1 )*SCALE
  244.          BTMP( 2 ) = BTMP( 2 )*SCALE
  245.       END IF
  246.       X2( 2 ) = BTMP( 2 ) / U22
  247.       X2( 1 ) = BTMP( 1 ) / U11 - ( U12 / U11 )*X2( 2 )
  248.       IF( XSWAP ) THEN
  249.          TEMP = X2( 2 )
  250.          X2( 2 ) = X2( 1 )
  251.          X2( 1 ) = TEMP
  252.       END IF
  253.       X( 1, 1 ) = X2( 1 )
  254.       IF( N1.EQ.1 ) THEN
  255.          X( 1, 2 ) = X2( 2 )
  256.          XNORM = ABS( X( 1, 1 ) ) + ABS( X( 1, 2 ) )
  257.       ELSE
  258.          X( 2, 1 ) = X2( 2 )
  259.          XNORM = MAX( ABS( X( 1, 1 ) ), ABS( X( 2, 1 ) ) )
  260.       END IF
  261.       RETURN
  262. *
  263. *     2 by 2:
  264. *     op[TL11 TL12]*[X11 X12] +ISGN* [X11 X12]*op[TR11 TR12] = [B11 B12]
  265. *       [TL21 TL22] [X21 X22]        [X21 X22]   [TR21 TR22]   [B21 B22]
  266. *
  267. *     Solve equivalent 4 by 4 system using complete pivoting.
  268. *     Set pivots less than SMIN to SMIN.
  269. *
  270.    50 CONTINUE
  271.       SMIN = MAX( ABS( TR( 1, 1 ) ), ABS( TR( 1, 2 ) ),
  272.      $       ABS( TR( 2, 1 ) ), ABS( TR( 2, 2 ) ) )
  273.       SMIN = MAX( SMIN, ABS( TL( 1, 1 ) ), ABS( TL( 1, 2 ) ),
  274.      $       ABS( TL( 2, 1 ) ), ABS( TL( 2, 2 ) ) )
  275.       SMIN = MAX( EPS*SMIN, SMLNUM )
  276.       BTMP( 1 ) = ZERO
  277.       CALL DCOPY( 16, BTMP, 0, T16, 1 )
  278.       T16( 1, 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
  279.       T16( 2, 2 ) = TL( 2, 2 ) + SGN*TR( 1, 1 )
  280.       T16( 3, 3 ) = TL( 1, 1 ) + SGN*TR( 2, 2 )
  281.       T16( 4, 4 ) = TL( 2, 2 ) + SGN*TR( 2, 2 )
  282.       IF( LTRANL ) THEN
  283.          T16( 1, 2 ) = TL( 2, 1 )
  284.          T16( 2, 1 ) = TL( 1, 2 )
  285.          T16( 3, 4 ) = TL( 2, 1 )
  286.          T16( 4, 3 ) = TL( 1, 2 )
  287.       ELSE
  288.          T16( 1, 2 ) = TL( 1, 2 )
  289.          T16( 2, 1 ) = TL( 2, 1 )
  290.          T16( 3, 4 ) = TL( 1, 2 )
  291.          T16( 4, 3 ) = TL( 2, 1 )
  292.       END IF
  293.       IF( LTRANR ) THEN
  294.          T16( 1, 3 ) = SGN*TR( 1, 2 )
  295.          T16( 2, 4 ) = SGN*TR( 1, 2 )
  296.          T16( 3, 1 ) = SGN*TR( 2, 1 )
  297.          T16( 4, 2 ) = SGN*TR( 2, 1 )
  298.       ELSE
  299.          T16( 1, 3 ) = SGN*TR( 2, 1 )
  300.          T16( 2, 4 ) = SGN*TR( 2, 1 )
  301.          T16( 3, 1 ) = SGN*TR( 1, 2 )
  302.          T16( 4, 2 ) = SGN*TR( 1, 2 )
  303.       END IF
  304.       BTMP( 1 ) = B( 1, 1 )
  305.       BTMP( 2 ) = B( 2, 1 )
  306.       BTMP( 3 ) = B( 1, 2 )
  307.       BTMP( 4 ) = B( 2, 2 )
  308. *
  309. *     Perform elimination
  310. *
  311.       DO 100 I = 1, 3
  312.          XMAX = ZERO
  313.          DO 70 IP = I, 4
  314.             DO 60 JP = I, 4
  315.                IF( ABS( T16( IP, JP ) ).GE.XMAX ) THEN
  316.                   XMAX = ABS( T16( IP, JP ) )
  317.                   IPSV = IP
  318.                   JPSV = JP
  319.                END IF
  320.    60       CONTINUE
  321.    70    CONTINUE
  322.          IF( IPSV.NE.I ) THEN
  323.             CALL DSWAP( 4, T16( IPSV, 1 ), 4, T16( I, 1 ), 4 )
  324.             TEMP = BTMP( I )
  325.             BTMP( I ) = BTMP( IPSV )
  326.             BTMP( IPSV ) = TEMP
  327.          END IF
  328.          IF( JPSV.NE.I )
  329.      $      CALL DSWAP( 4, T16( 1, JPSV ), 1, T16( 1, I ), 1 )
  330.          JPIV( I ) = JPSV
  331.          IF( ABS( T16( I, I ) ).LT.SMIN ) THEN
  332.             INFO = 1
  333.             T16( I, I ) = SMIN
  334.          END IF
  335.          DO 90 J = I + 1, 4
  336.             T16( J, I ) = T16( J, I ) / T16( I, I )
  337.             BTMP( J ) = BTMP( J ) - T16( J, I )*BTMP( I )
  338.             DO 80 K = I + 1, 4
  339.                T16( J, K ) = T16( J, K ) - T16( J, I )*T16( I, K )
  340.    80       CONTINUE
  341.    90    CONTINUE
  342.   100 CONTINUE
  343.       IF( ABS( T16( 4, 4 ) ).LT.SMIN )
  344.      $   T16( 4, 4 ) = SMIN
  345.       SCALE = ONE
  346.       IF( ( EIGHT*SMLNUM )*ABS( BTMP( 1 ) ).GT.ABS( T16( 1, 1 ) ) .OR.
  347.      $    ( EIGHT*SMLNUM )*ABS( BTMP( 2 ) ).GT.ABS( T16( 2, 2 ) ) .OR.
  348.      $    ( EIGHT*SMLNUM )*ABS( BTMP( 3 ) ).GT.ABS( T16( 3, 3 ) ) .OR.
  349.      $    ( EIGHT*SMLNUM )*ABS( BTMP( 4 ) ).GT.ABS( T16( 4, 4 ) ) ) THEN
  350.          SCALE = ( ONE / EIGHT ) / MAX( ABS( BTMP( 1 ) ),
  351.      $           ABS( BTMP( 2 ) ), ABS( BTMP( 3 ) ), ABS( BTMP( 4 ) ) )
  352.          BTMP( 1 ) = BTMP( 1 )*SCALE
  353.          BTMP( 2 ) = BTMP( 2 )*SCALE
  354.          BTMP( 3 ) = BTMP( 3 )*SCALE
  355.          BTMP( 4 ) = BTMP( 4 )*SCALE
  356.       END IF
  357.       DO 120 I = 1, 4
  358.          K = 5 - I
  359.          TEMP = ONE / T16( K, K )
  360.          TMP( K ) = BTMP( K )*TEMP
  361.          DO 110 J = K + 1, 4
  362.             TMP( K ) = TMP( K ) - ( TEMP*T16( K, J ) )*TMP( J )
  363.   110    CONTINUE
  364.   120 CONTINUE
  365.       DO 130 I = 1, 3
  366.          IF( JPIV( 4-I ).NE.4-I ) THEN
  367.             TEMP = TMP( 4-I )
  368.             TMP( 4-I ) = TMP( JPIV( 4-I ) )
  369.             TMP( JPIV( 4-I ) ) = TEMP
  370.          END IF
  371.   130 CONTINUE
  372.       X( 1, 1 ) = TMP( 1 )
  373.       X( 2, 1 ) = TMP( 2 )
  374.       X( 1, 2 ) = TMP( 3 )
  375.       X( 2, 2 ) = TMP( 4 )
  376.       XNORM = MAX( ABS( TMP( 1 ) )+ABS( TMP( 3 ) ),
  377.      $        ABS( TMP( 2 ) )+ABS( TMP( 4 ) ) )
  378.       RETURN
  379. *
  380. *     End of DLASY2
  381. *
  382.       END
  383.